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Abstract 

The Schrodinger equation (Hip)(r) = {E + u^W {r))if}{r) with an energy-dependent complex 
absorbing potential — ueW(t), associated with a scattering system, can be reduced for a special 
choice of ue to a harmonic inversion problem of a discrete pseudo-time correlation function y(t) = 
(j^U 1 ^. An efficient formula for Green's function matrix elements is also derived. Since the exact 
propagation up to time 2t can be done with only ~ t real matrix-vector products, this gives an 
unprecedently efficient scheme for accurate calculations of quantum spectra for possibly very large 
systems. 
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Complex absorbing potentials. The spectrum of a quantum scattering system can 
be characterized by solving the boundary value problem associated with the Schrodinger 
equation, 

{H^){r) = Eip(r). (1) 

The bound states then have real energies E, with solutions ip(r) exponentially localized in 
space. The resonance states (Siegert states |JJ) have complex energies with ImE < 0. They 
behave like bound states in some compact subset Q of the configuration space, but eventually 
grow exponentially outside of Q, due to the outgoing asymptotic boundary conditions. 

By introduction of a so-called optical (or absorbing) potential — uW(r) with Imu > 
and real W(r) > 0, that vanishes for r G Q and smoothly grows outside fl, the solutions 
ip(r) of 

(Hip)(r) = (E + uW(r))ip(r) (2) 

are damped outside Q, the physically relevant region ||. This forces them to behave like 
bound states everywhere without significantly affecting the energies E. In this framework 
the physically relevant part of the system is, therefore, dissipative and satisfies (0) only 
for r G fl Moreover, a general multichannel scattering problem can be considered with a 
numerically convenient form of W(r), independent of the choice of coordinate system. 

Although to satisfy ImE < one only needs Imw > 0, traditionally one simply puts 
u = i and gets the nonhermitian eigenvalue problem (H — iW)tf) = Eip. The latter is 
generally much easier to handle numerically than the boundary value problem ([!]). As we 
shall see, energy-dependent choices u = ue are particularly useful. 

The introduction of the absorbing potential leads to the damped Green's function |3|-f| 

G W {E) := (E-H + UeW)- 1 . (3) 

Under suitable conditions on UEW(r), it is probably possible to prove, similar as in ref. 
for the traditional case ue = i, that G\y{E) converges for any real E (and also for ImE > 0) 
weakly to the ordinary Green's function 

G(E) = \im(E -H + is)' 1 . 

Practically, one usually needs to evaluate only certain matrix elements (f T G{E)%p ) the basic 
numerical objects of quantum physics from which most other quantities of interest (scattering 
amplitudes, reaction rates, etc.) can be computed (see, e.g., refs. 0,f|,|6|,[7j ) . If both and ijj 
have support in Q, they are well approximated by <j) T Gw{E)ip. 

Unfortunately, for very large systems with high density of states one may encounter 
numerical difficulties when trying to diagonalize a large nonhermitian matrix H — ueW 
or solve the linear system (E — H + ueW)X(E) — ip at many values of E using general 
iterative techniques for nonhermitian matrices. However, one can devise alternative iterative 
techniques to solve (|]) by exploiting the special structure of the quantum scattering problem. 

From now on, we assume that the Hilbert space is discretized so that the states are 
vectors ip G C K and if, W are real symmetric K x K matrices, W diagonal, as, e.g., in the 
case of a discrete variable representation M. 
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To simplify the following equations we further assume without loss of generality that 
the discretized Hamiltonian is shifted and scaled so that < 1 for any state ip, where 

we have defined the expectation value as (H)^ := ip*Hip /%p*%p and * denotes conjugate 
transposition. Such scaling is implemented routinely in the framework of the Chebyshev 
polynomial expansion. 

We now consider the special choice 

1+1,2 

u E = E + iVT^E 2 ~, or%) = ^-, (4) 

which after insertion into @ gives a nonlinear eigenvalue problem for u = ue, more useful 
than (||): 

1 + u 2 D 

Htp = — y- — ip with D = 1 + 2W. (5) 

We may think of this equation as an eigenvalue problem 

Hip = E{u)ip (6) 
involving an operator-valued li-dependent energy 

„, . 1 + u 2 D , . 

E(u) = (7) 

that reduces in Q (where D(r) = 1) to the constant (f|). 

In ref. a similar nonlinear eigenvalue problem with E(u) = (D~ x + u 2 D)/2u was 
implicitly encountered leading to a numerical scheme to compute, e.g., the complex resonance 
energies by harmonic inversion of a "discrete-time" correlation function, generated by a 
damped Chebyshev recursion. Here, related results are derived rigorously and in a more 
general framework. In particular, a pseudo-time Schrddinger equation is derived that allows 
one to achieve a substantial numerical saving compared to the previous works. 

The eigenpairs (u k , ip k ) of @ can be used to evaluate the physically interesting quantities 
(e.g., the complex resonance energies E k , scattering amplitudes, etc.). However, because of 
the nonlinearity they have somewhat different properties from those of the regular nonher- 
mitian eigenvalue problem, which we now proceed to derive. 

Theorem 1: Completeness and Orthogonality. The nonlinear eigenvalue problem @ 
has at most 2K distinct eigenvalues u k - Thus, there are at most K physical eigenvalues E k 
of (D with \m.E k > 0. 

If there are 2K distinct eigenvalues u%, . . . ,U2K with associated eigenvectors tp k satisfying 
(H), then for any vector (j) a G C K there is a set of 2K numbers 9 ak satisfying the completeness 
relations 

2K 2K 
0a = °ak1pk, = d ak U k 1p k . (8) 
k=l k=l 

Furthermore, ip k satisfy the orthogonality relations 
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■0j (1 - UjUkD)^ = 5 jk (9) 

for all j ^ k. If the eigenvectors can be normalized such that (|9p also holds for j = k, then 
(|) holds with 

Oak = iflK (10) 

Proof. Let 
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7 N 

-D- 1 D- l {2H) 



where 7 denotes the K x K unit matrix. The ordinary eigenvalue problem 

U$ = vi (11) 

yields, with ip = Qj) , 

ip' = uip, D' 1 ^ + u 2 ij -uD-\2H)ij = 0. 
After multiplication by D/2u, we find (d). In particular, the eigenpairs (uk,^) °f ^7 satisfy 

4 = C ^ fc ^ (12) 

where (u k ,ip k ) is an eigenpair of ([5|). Since conversely, any such eigenpair determines an 
eigenpair of U, the nonlinear eigenvalue problem fl5|) has at most distinct eigenvalues. If 
there are 2K distinct eigenvalues ui, . . . , u 2 k, the matrix U is diagonalizable, and there is a 
basis ipi, ... , ip2K of eigenvectors of the form (0). Therefore, we may write 



with uniquely determined coefficients 9 a k- This gives (^). 

Using (H) and the symmetry of H and 7? we may compute ipjHi(j k in two different ways: 



lT l+u 2 k D T l+u 2 D 



For j 7^ /c, by assumption Uj ^ Uk- Thus, we may take difference, multiply by 2ujUk/ (uk—Uj), 
and find (0). For j = fc, we may achieve by normalizing the eigenvectors, provided that 
the left hand side of (§) does not vanish. Using @ and (§) we find (|10[) : 



Oak = J2j®<xi S jk = JLjQaj^Ji 1 ~ UjU k D)tlj k 
= &k = Ipk&a- 
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We now consider an eigenpair (u, ip) of fl5|) with ^ 0. Multiplying (j5|) by 2uip* gives 
the quadratic equation 



u {D)^ — 2u(H)^ + 1 = 0. (13) 



The solutions of fll3|) are 



»= k ■ (14) 

Since (D) = 1 + 2(W)^ > 1 and \ {H)^\ < 1, the square root is real and 

■■j2 m+(Du-(H)i i 

H " w "w;- 1 - (15) 

Thus, u is a complex number lying in the unit disk. Moreover, \u\ = 1 iff (W)^ = 0, i.e., iff ip 
has support in fl, which is the case for the bound states. The states with \u\ ~ 1 correspond 
to the narrow resonances. Due to ([14]) the solutions of (|5p come in complex conjugate pairs 
(u,ip) and (u,ip). The physically relevant eigenenergies with ImE < come from u with 
Imu > 0. 



Note that a similar analysis of a quadratic eigenvalue problem was carried out in ref. [JTU 



arising from the use of the Bloch operator L, rather than an absorbing potential. There, 
equation @ is considered with E(u) = [iuL + u 2 I)/2, and u is a momentum variable close 
to the real axis, instead of a number close to the unit circle. However, this equation can 
only be used for less general, single- channel scattering problems and, besides, it is hard to 
solve efficiently using iterative techniques. 

Reduction to a harmonic inversion problem. Consider the pseudo-time Schrddinger 
equation defined by the recurrence 

cj)(t) =D- 1 (2H(f)(t-l)-(f)(t-2)) (t = 2,3,...) (16) 

with 0(0) = 0o and (f)(1) = 0. (This choice of initial conditions is most convenient, although 
other initial conditions with 0(1) ^ yield analogous results. A similar 3-term-recurrence 
with another choice of special initial conditions leading to "modified Chebyshev recurrence" 
was considered in refs. PJTTf.) Since D is diagonal and matrix- vector products H<p are 
usually cheap to form, 0(0), . . . , 0(T) are computable using O(KT) operations and a few 
vectors stored at a time. If the initial vector 0o is real, only real arithmetic is needed. 
By Theorem 1, we can write (|i"6|) as 

and, therefore, 

2K 

0W = E^4^- (17) 





k=l 
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This power expansion is very important, and is analogous to the (physical time) expan- 
sion (f)(t) = e~ ltH (fi(0) = Y,k=i@ke~ ltEk, 4>k for the solutions of the standard time- dependent 
Schrodinger equation. It allows one to reap all the benefits of time-dependent methods (see, 



e.g., refs. [p!2| , |13| , p|) without having to deal with the time-dependent Schrodinger equation, 
which is hard to solve accurately at long times t in the case of nonhermitian Hamiltonian. 
Instead, only the much more benign and numerically very stable equation (|16|) must be 
solved. 

By ([171), the pseudo-time cross-correlation function 

y a {t) := fy{t) (t = 0,l...) (18) 

of a state <ft a has the form 

2K 

y a (t) = dakul (19) 



k=l 



with 



d ak = (j^lpk^k^a = ok 9 ak . (20) 



This reduces the nonlinear eigenvalue problem (|5[) to solving the harmonic inversion problem, 
i.e., to finding the spectral parameters (uk,d a k) {k = l,...,2K) satisfying ([19|) for the 
sequence y a (t) computed by (|T6| ) and ([18]). Since by (|15|) the sequence y a if) is bounded 
and the spectral mapping (J|) moves the physically relevant eigenvalues u k close to the unit 
circle, this is an efficiently tractable problem, even in very large dimensions p!^j9||. 



Time doubling of an autocorrelation function. As is well known, a true time auto- 
correlation function at time t can be computed by solving the time-dependent Schrodinger 
equation up to time t/2, since one can use 

C(t) := T e"^V = (e-^ 2 0) T (e-^/ 2 0). 



For the Chebyshev autocorrelation function c(t) := 0„0(t), based on ( |T6D with D = I and 



the initial conditions 0(0) = 0o> 0(1) = H(j) , a factor of two saving is also well known (see, 
e.g., the discussion in ref. 0): 

c{2t) = 20(t) T 0(t) - c(0), c{2t + 1) = 20(t) T 0(t + 1) - c(l). 



This expression was used in |T5[ for resonance computation implementing a damped Cheby- 
shev recursion. However, being approximate, it only worked for sufficiently narrow reso- 
nances. In the present framework, we can write the pseudo-time cross-correlation function 

as 

Va{t) = 

which suggests that an exact doubling scheme exists. This is now derived for the autocor- 
relation function as only O is propagated. 
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Theorem 2: The Doubling Scheme. For vectors <fi(t) and 0(s) satisfying the pseudo- 
time Schrodinger equation (|16|) with initial conditions 0(0) = 0o, 0(1) — 0, the autocorrela- 
tion function yo(t) := 0q 0(t) satisfies 

y (s + t) = 0(s) T 0(t) - 0(s + l) T D0(t + 1) =: z(s, t). 
Proof: This follows from the power expansion ( ]T7| ) and the orthogonality relations ([5|): 

2K 

z(s,t) = ^ OojOokUjuKiffijjk - UjUkijjD^k) 
j,k=i 

2K 

= Y,vZkU s k +t = yo(s + t). 
k=l 

Hardly any additional storage will be needed if the sequence yo(t) (t = 0, 2T — 2) is 
generated by 

y (2t) = z(t,t), y (2t-l) = z{t,t-l), (21) 

concurrently with the computation of 0(t) using t = 0, ...,T. In exact arithmetic the har- 
monic inversion of the doubled sequence yo(t) will give the exact results if T > 2K, thus, 
using only T ~ IK of matrix-vector products. However, this is impractical as it would 
formally require to solve a T x T eigenvalue problem. To reduce the computational burden 
and to maintain numerical stability the eigenvalues are extracted very efficiently in a small 
Fourier subspace by the Filter Diagonalization Method P^fl. In this case, the required 



length 2T of the doubled sequence needed to converge an eigenenergy E k (cf. Eq. |) will be 
defined by the locally averaged density of states p(E) for E k ~ E 0. 

Theorem 3: The Green's function matrix elements. Under the assumptions of The- 
orem 1, let 0(t) be a solution of the pseudo-time Schrodinger equation (|l^) with initial 
conditions 0(0) = 0o, 0(1) = 0. Then the matrix elements of the damped Green's function 
(|) with u E = E + — E 2 are 



2K 



dj3 k d ak 2uEU k 



Gw(E)<p a = E T , (22) 



k=l 



d ok u k - u E 



where the three sets of spectral parameters {dp k , u k }, {d ak , u k } and {d ok , u k } (with identical 
eigenvalues u k ) satisfy the harmonic inversion problem (|19D for the cross-correlation functions 
yp{t), y a {t) and y (t), respectively. 

Proof. For an eigenpair (u k ,ip k ) of (^|) we can write 

(E-H + u E W)^j k = u *~ UE (i - Uk u E D)ip k . 

2u k u E 

Then (|) and ([TUp imply 
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(E-H + u E W) \ 

^ u k -u E 

= Qak^Pk ~ U E D I ^ 9 ak U k 4) k ) = ^ 0Q, fc Vfc = </><*• 
fc V k Ik 

Multiplying this by (f>J(E — H + ueW)' 1 and using 6p k = 4>pip k we obtain 

4 > ljGw{E)(j) a = dgkOgk , (23) 

^ u k - U E 

Now replacement of dp k d ak by dp k d ak /do k gives (^). 

Note that (^) also gives an explicit expression for the damped Green's function in terms 
of the eigenpairs: 

G w (E) = f:^^^l (24) 
k% u k - u E 



A formula similar to (|22|) was obtained (without a rigorous derivation) in ref. |16 



in 



the framework of the damped Chebyshev recursion in place of (|16|). However, here, due to 
the doubling scheme fl2T|), only half the number of matrix- vector products will be needed to 
obtain the same amount of information. 

Thus, we have a very efficient and stable method to extract the complete spectral and dy- 
namical information of a general (multichannel) quantum scattering system using a minimal 
number of matrix- vector products. This will be demonstrated numerically in a forthcoming 
publication. 
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